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1.  Introduction 


fn  this  paper  we  present  the  Modified  ConjuRate  Residual  (MCR) 

^ * 

Method,  a stabilized  version  of  Luenberger's  Method  of  Conjugate 
Residuals  , for  solving  large  sparse  systems  of  linear  equations. 

This  Iterative  method  has  special  significance  when  the  system  Is  not 
positive  definite  so  that  methods  like  Conjugate  Gradients  (2]  are 
inapplicable.  In  the  special  case  when  the  system  is  positive  definite, 
MCR  reduces  to  one  of  the  family  of  general  conjugate  gradient  methods 
discussed  by  Hestenes^  {^] . In  section  2 we  present  the  MCR  method  and 
its  properties.  In  section  S we  derive  general  error  bounds,  and  In 
section  4 we  present  an  efficient  computational  version  of  the  MCR 
algorithm.  In  section  5 we  Introduce  a model  Indefinite  problem,  simple 
finite  difference  approximations  to  the  forced  vibration  problem  In  two 
and  three  dimensions  and  apply  the  convergence  results  of  section  3 to 
obtain  error  bounds.  Although  these  problems  are  very  specialized, 
nearly  all  our  results  are  applicable  to  the  solution  of  the  linear 
systems  arising  from  the  application  of  finite  difference  or  finite 
element  methods  to  more  general  self-adjoint  second  order  elliptic 
partial  differential  equations  In  more  general  domains.  Finally  in 
section  6 we  describe  the  results  of  numerical  experiments  on  the  model 
problem.  Proofs  of  the  theoretical  results  and  further  details  will 
appear  in  [1]. 


□ C 
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2.  The  Hodlfled  Conlugate  Residual  (MCR)  Method 

Consider  the  NxN  linear  system 

(2.1)  Ax  - f 

where  A Is  non-singular  and  symmetric.  To  approximate  x.  the  MCR  method 

generates  a sequence  of  direction  vectors  p^  which  are  mutually 
2 

A -orthogonal,  and  computes  the  Iterates  x^.Xi.x^,...  by  moving  at 
8 L 

the  ( l-fl ) step  In  the  direction  of  p^  so  that 

Xj^^  «■  Xj^  + where  Is  chosen  to  minimize  the  error 

functional 

(2.2)  E(a^)  - (f-Ax^^^,  f-Ax^^^). 


Due  to  the  special  way  In  which  the  direction  vectors  are  chosen  the 
unl-dlrectlonal  minimizations  correspond  to  minimization  In  the  whole 
subspace  spanned  by  the  p^'s,  so  that  the  obtained  actually 

minimizes  E(x)  on  the  affine  subspace  x^  + (Pq.P j,P2» • • • »Pj) • 
where  {...}  denotes  the  subspace  spanned  by  the  vectors  enclosed  In 
the  brackets. 


In  particular,  we  use  the  following  method,  based  on  the  Lanczos 
algorithm  [7]  for  generating  the  direction  vectors: 


Choose  Pq,  define  p ^ 0,  and  for  1~0,1,...,  let 


(2.3a) 


APi  - - «jPi_i. 


where 

(2.3b) 


2 

- (A  p^,Ap^)/(Ap^,Ap^)  and 

«1  - (A^P4*APj_j)/Ap^_j,Ap^_j). 
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We  can  show  that 

(2.4,  (Pj.aSj)  - 0.  1 O. 

We  now  consider  the  following  scheme  for  computing  the  approximations 
Xj  to  x: 


Given  an  arbitrary  Initial  guess  x^  to  x,  compute 
Cq  " f-AXp,  set  pQ  " Cq  and  for  1 - 0,1,2,...  compute 


(2.5) 


®1  “ ('^i*APj)/(APj.Ap^) 


’‘l.l  ■ "l  ^ “l"! 

■■l.l  ■ "l  - “1*^1 

and  generate  p^^^  by  (2.3), 


In  section  3 we  show  that  the  above  procedure  can  be  made  more  efficient 
in  certain  cases.  However,  for  now,  we  refer  to  (2.5)  as  the  MCR 
algorithm  and  proceed  to  outline  Its  properties  and  convergence  results. 


We  can  show  that  the  following  relations  hold  for  the  MCR  method  of 


(2.5): 

, 

(2.6a) 

APj  C (pQ»Pj^»  • • • »Pj^J  ^ 

(2.6b) 

*^1  " <VPi Pi^ 

(2.6c) 

{Po*Pi' •••»Pi>  ■ ^Pq’^O’ 

(2.6d) 

(r^.Apj)  - 0,  j < 1, 

(2.6e) 

(r^.ACj)  - 0,  1 y j. 

(2.6f) 

(fj.Ap^)  - (r^,Ap^),  j < 

(2.6g) 

and  r^  y C “>  p^  y 0 

t 

(0  denotes  the  zero 

(ro.Arg. 


.A^o> 
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Ue  also  have  the  result: 

Theorem  2. 1;  For  each  1,  minimizes  E(x)  over  the  affine 

subspace  + (Pq.Pj Pj>* 
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3.  Error  Bounds 

The  fundamental  minimization  property  of  Theorem  2.1  allows  us  to 
view  the  MCR  method  as  one  In  which  the  error  functional  Is  minimized 
over  subspaces  of  increasing  dimension,  and  thereby  to  prove  that  It 
converges  In  at  most  N steps  to  the  unique  solution  of  Ax  ■ f. 

Moreover,  we  can  obtain  error  bounds  by  a procedure  analogous  to  the  one 
used  for  Conjugate  Gradients  (3). 

Using  Theorem  2.1  and  property  (2.6c)  «/e  can  write 
^ j 1 

X,  , “ X-  + I s.A  r where  {s,}.  are  scalars 
i+1  0 j,0  J ° 

■ Xq  + P^(A)  r^  where  P^(A)  is  a polynomial  of  degree  at  most  1 in  A. 

Then, 


x - x^^j  - (I  - P^(A)A)  (X  - Xjj) 

and 

(3.1)  E(x  - ) ■ (f  - f 

(Rl^l(A)ro,Ri_^l(A)ro) 

where 

Ri^l(A)  - I - P^(A)A. 

Let  us  define  R^  to  be  the  set  of  polynomials  of  degree  at  most 

m in  y satisfying  R (0)  “ 1.  Then,  due  to  (3.1)  and  Theorem  2.1,  we 
n 

sL 

can  visualize  MCR  as  a method  that  chooses,  at  the  (1+1) 

iteration,  the  particular  polynomial  R e R that  minimizes  the 

m m 

functional  E(x) , l.e.. 


E(Xi^l) 


min 
R 


1+1  \+l 


i 


i 
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k. 


Since  A Is  symmetric,  there  exist  orthonormal  eigenvectors  v ^ , 
J ••  of  A such  that 


AVj  ■ ^ " 1*2, 


where  {X.}.  , are  the  eigenvalues  of  A. 

J J-1  N 

” N 

We  have  r^  - £ t.v.  for  some  scalars  (t.>.  .. 

0 j.i  J J j J-1 


Then 


N 


N 


and  N N 

"iti  ' Vi  ^ ‘ J ' 

^22 

- min  £ t,  R.  , (X  ) 

R...  c J-1  J J 


i+1  1+1 


>[ 


min 

^+1  ' ^1+1 


max  lRi+i(>j)l|^ 


(3.2) 


Vl 


where  Q - min  max 

"ixl  ‘ ^+1  ^ 


Depending  upon  the  spectrum  of  A and  what  is  known  about  it,  we  can 
use  different  techniques  to  get  bounds  on  • W®  consider  two 
approaches.  First,  we  treat  the  case  in  which  at  roost  a few  eigenvalues 
lie  on  one  side  of  the  origin  and  we  know  the  end-points  of  the  interval 
containing  the  eigenvalues  on  the  other,  for  example,  in  the  solution  of 
difference  approximations  to  the  forced  vibration  problem.  For  this 
case  we  have  the  following  theorem: 


Theorem  3j_l:  If  the  eigenvalues  of  the  matrix  A lie  in 
{X,,X  ,...,X  } u [a,b]  where  1 £ m < N,  a,  b>0,  and  a < 0 


-7- 


for  i £ i ^ then  for  k ^ 0 the  Iterates  x given  by  MCR 
^ It 

satisfy 

(3.3)  II  t - II  j < II  « - A.,  II  j 

and 

(3.4)  II  X - II2  < 2C  /7(K^)  II  X - Xq  II  2 

m ^ -b 

where  a - a/b,  C ■ n (— J — ) is  a constant  Independent 

2 J-1  J 2 

of  k,  and  x(a  ) denotes  the  condition  number  of  A . 

If  the  matrix  A Is  positive  definite,  (3.3)  and  (3.4)  hold  with 
C - 1,  and  we  have  the  same  rate  of  convergence  for  Modified  Conjugate 
Residuals  as  for  Conjugate  Gradients  (cf.  [3]). 


In  the  second  approach  we  assume  that  all  we  know  is  the  end-points 
of  the  two  intervals,  one  on  either  side  of  the  origin,  that  contain  the 
eigenvalues  of  A.  Using  Lebedev's  results[S]  we  can  show  the  follo%rlng 
result. 


Theorem  3. 2;  If  the  eigenvalues  of  the  matrix  A lie  in  (-a,-b)  (c,d) 

where  a,b,c,d  > 0 and  a-b"d-c>0  then  for  i>0  the  iterates 
given  by  MCR  satisfy 


i+1 


and 


f - Ax 


1+1  "2 


U±1 

. ..,1  - /I.  L 2 J,.  , 

<2(7-^)  Ilf-Axp 

1¥J 


A - Ai.i  "2^2 


X - X 


0 "2 


2 2 
where  B • (be) /(ad)  and  ic(A  ) denotes  the  condition  number  of  A • 


j 
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We  note  that  Theorem  3.2  requires  that  the  two  Intervals  on  either 
side  of  the  origin,  know  to  contain  the  eigenvalues  of  A be  of  equal 
length.  If  we  have  that  a-b^d-cwe  can  take  the  smaller  of  the 
two  Intervals  and  extend  It  away  from  the  origin  to  get  Intervals  of 
equal  length  on  either  side  of  the  origin  and  then  apply  Theorem  3.2  to 
obtain  bounds  on  the  rate  of  convergence. 
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4.  The  MCR  Alaorlthm  - Cowputatlonal  Version 

In  this  section  we  show  that  in  certain  cases  we  can  reduce  the 
cost  of  an  iteration  of  the  HCR  method  by  using  a less  expensive  process 
than  (2.3)  for  generating  the  new  direction  vector.  We  present  the 
computational  version  of  the  MCR  method  In  Its  reduced  form  for  both 
positive  definite  A and  general  symmetric  A.  The  algorithms  are 
numerically  stable,  easy  to  program  and  require  very  little  storage  and 
work  per  iteration. 


In  particular,  the  following  identity  holds. 


Theorem  4. 1 ; If  a^  i*  0, 


(4.1.)  fl.l  - ■'l.l  ♦ Vl- 

where 

(4.1b)  b^  - (-Ar^^j,Ap^)/(Apj,Ap^). 

and  Is  given  by  relations  (2.3),  then 


i+1 


-^Pl+1 


The  above  theorem  shows  that  if  a^  ^ 0,  equations  (4.1)  generate 
the  same  direction  as  the  more  expensive  (2.3).  Only  the  normalization 
is  different,  but  since  that  does  not  effect  the  validity  of 
orthogonality  relations,  it  is  clear  that  all  the  results  of  sections  2 
and  3 hold. 


In  the  special  case  of  positive  definite  A,  we  can  show  (cf.  (!]) 
that  a^  is  never  zero  and  that  (4.1)  can  be  used  to  generate  the 
direction  vectors  at  every  Iteration.  We  also  use  alternative  formulae 
for  a^  and  b^  summarized  in  (4.2): 


A 
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1- 


MCR  algorithm  positive  definite) : 


Choose  Xq,  compute  ■ f - Ax^,  set  Pq  ■ and  for  1 
0, 1 , 2, . . > compute 


(4.2) 


®1  " 


*1+1  “ *1  ^ ^Pl 

"l+l  “ "l  - "I'^l 


1+1 


‘■l+l  ^ ‘’iPl 


The  number  of  multiplications  per  Iteration  for  the  above  procedure  Is 


5N  + 2 plus  a matrix-vector  product,  viz.  Ar^^^;  Ap^^^j  Is  computed  as 
*^1+1  ■ “wi  * 


Storage  requirements  are  5N  plus  whatever  Is  required  to  store  the  upper 
triangular  part  of  A. 


For  general  symmetric  A,  we  could  have  a^  equal  to  zero  In  which 
case  the  longer  Iteration  (2.3)  would  have  to  be  used.  However,  each 
time  that  a^  0 we  could  use  the  shorter  Iteration  and  save  some  work. 
In  practice,  though,  due  to  finite  precision  and  roundoff  limitations, 
it  may  not  be  possible  to  decide  whether  a^^  Is  zero  or  not.  We  can 
show  that  If  a^  “ 0 and  we  use  the  shorter  version  (4.1)  to  generate 
the  direction  vector,  the  algorithm  may  get  tuck  at  a point  and  serious 
instabilities  may  occur.  Therefore,  we  decide  the  use  the  shorter 
version  only  when  la^^l  > e where  e > 0 is  some  threshold  chosen  a 
priori  and  large  enough  such  that  roundoff  and  precision  will  not  cause 
us  to  choose  (4.1)  when  a^  « 0.  We  note,  however,  that  if  0 < a^  < c 
we  are  using  the  longer  Iteration  at  step  1,  even  though  the  shorter  one 


L 


Is  sufficient. 


We  can  also  show  that. 


- c^{Ap^,Ap^)/(Ap^_^,Ap^_j ) 


where 


1 if  I a^_j I < c 
-l/ai_^  if  I > £ 


We  thus  have  the  following  computational  version  of  the  Modified 
Conjugate  Residual  Method  for  any  symmetric  indefinite  linear  system 


Ax  = f: 


MCR  algorithm  (A  symmetric) : 


Choose  Xq,  and  e > 0;  compute  r^  - f - Ax^^,  set  p^  « r^  and 
for  i =•  0,1,2,  ...  compute 


Hf  - (r^,Ap^)/(Ap^,Ap^) 


*1+1  ■ *1  “l"! 


*1+1  ■ *1  - ■>! 


If  |a^|  < e , 


(4.3) 


Pi+1  * ^Pi  - ^Pi  - «iPi-l 

2 

where  = (A  p^ ,Ap^) /(Ap^ ,Ap^) 
and  = c^(Ap^,Ap^)/(Ap^_^,Ap^_j) 


where  c . 


1 if  |aj_jl  < o 

-1/^.1  if  'Vl'  " " 


If  |a^|  > e. 


Pi+i  ■ *1+1  * *’i'’i 


where  ,Ap  )/(Ap^ , Ap  ) 
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Each  Iteration  of  (A. 3)  requires  only  one  matrix-vector  product. 

To  see  this,  suppose  that  we  have  already  computed  and  stored  the 

vectors  x^,  r^,  p^,  Ap^,  p^_j  and  Ap^_j  and  the  scalar 

(Ap^  1*^^!  1^*  ® only  matrix-vector  product 

St. 

performed  at  the  (l+l)  Iteration  is  A(Apj)  required  In  the 
computation  of  Ap^^^  which  Is  done  as  follows: 


^1+1  ■ - 'l^^l  “ *1^1-1 


If  |a^|  > w,  then  Ar^^^  Is  computed  and  Ap^^^  Is  obtained  as 


Ap 


l+l 


Ar 


l+l 


+ b^Ap^. 


It  Is  easy  to  see  that  In  addition  to  the  matrix-vector  product,  the 
Iteration  requires  7N  + 2 multiplications  If  [a^^l  > c and  9N  + 4 
multiplications  of  la^l  £ c. 


Storage  Is  required  for  the  six  vectors  x^,  r^,  p^,  Ap^,  p^  ^ 

and  Apj^  , each  of  which  gets  overwritten  by  Its  successor  at  the 

s t 2 

(1+1)  Iteration,  and  a seventh  vector  that  stores  A p^  If 

|a^|  < E and  Ar^^j^  If  |a^|  > c..  The  matrix  A Is  not  modified 

during  the  Iteration  process.  Since  It  Is  symmetric,  only  the  upper  (or 

lower)  triangle  may  be  stored.  If  A Is  sparse,  sparse  storage  schemes 

may  be  used. 


We  note  here  the  close  connection  between  MCR  and  Luenberger's 

method  of  Conjugate  Residuals  [6].  Indeed,  taking  e - 0 we  can  rewrite 

the  algorithm  (4.3)  to  give  Luenberger's  method.  In  this  case,  at 

Iteration  k + 1 we  are  required  to  answer  the  question  — "Is  a ’ 0?" 

k 

As  pointed  out  earlier,  this  decision  is  not  easy  In  the  presence  of 
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roundoff  error  so  that  Luenberger's  method  suffers  from  certain 
unresolved  computational  difficulties.  By  choosing  c sufficiently 
large  we  have  the  MCR  method  which  Is  free  from  these  problems. 
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5.  The  Model  Problems 

In  this  section  we  Introduce  model  problems  in  two  and  three 
dimensions — simple  finite  difference  approximations  to  Indefinite 
self-adjoint  second  order  elliptic  partial  differential  equations  In  a 
square  and  cube  respectively  (the  forced  vibration  problem). 

In  two  dimensions  we  consider  the  problem: 

(5.1a)  -Aw(x,y)  - ow(x,y)  - g(x,y),  (x,y)  e D 5 (0,1)  x (0,1) 

with  Dirichlet  boundary  conditions 
(5.1b)  w(x,y)  =•  0,  (x,y)  e aD 

where  u Is  a constant. 

I 

i 

I 

To  approximate  the  solution  to  this  problem,  we  cover  the  domain  D j 

with  a uniform  mesh  with  mesh-width  h = l/(n+l)  and  seek  a mesh  ' 

function  W(i,J)  which  is  an  approximation  to  w(lh,jh)  for  each  1 £ i, 
j £ n.  If  we  replace  the  differential  operator  by  the  familiar 
five-point  difference  approximation  at  each  interior  mesh  point, 
cf.  [8],  we  obtain  the  system  of  linear  equations 

(4-ah^)  W(1,J)  - W(i,J-l)  - W(l,j+1)  - W(i-l,j)  - W(l+l,j)  - h^G(l,J), 

(5.2)  1 < i.  J<  n, 

where  W(i,j)  » 0 if  i - 0 or  n+  1 or  if  j « 0 or  n+1  and  G(1,J)  » 
g(lh,jh). 


If  the  unknowns  W(l,j)  are  ordered  in  the  natural  row-by-row 
fashion,  the  system  (5.2)  is  reduced  to  a nxn  block  trldiagonal  system 
of  linear  equations 


I 
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where  T is  the  nxn  trldiagonal  matrix 


and  I Is  the  nxn  Identity  matrix.  The  matrix  A has  five  non-zero 
diagonals  and  is  symmetric. 

In  three  dimensions,  we  consider  the  problem: 

-Aw(x,y,z)  - aw(x,y,z)  - f{x,y,z),  (x,y,z)  e D = (0,1)  x (0,1)  x (0,1) 

with  Dirlchlet  boundary  conditions 

w(x,y,z)  » 0,  (x,y,z)  e oD 

where  o is  a constant. 
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B -I  0 


where  B is  a n xn  matrix  of  the  form 


D is  the  nxn  tridiagonal  matrix. 


h « l/(n+l)  and  the  I's  denote  identity  matrices  of  suitable  order. 

If  a £ 0,  the  corresponding  matrices  A for  the  two  and  three 

dimensional  model  problems  are  positive  definite.  However,  if  a is 

sufficiently  positive,  the  corresponding  matrices  are  strictly 

indefinite,  i.e.,  they  have  both  positive  and  negative  eigenvalues.  We 

2 

also  have,  in  the  notation  of  Theorem  3.1,  a = 0(h  ).  Applying 
Theorem  3.1,  we  obtain  the  following  convergence  results: 

Theorem  5.1:  For  the  model  problems,  MCR  requires  0(n  log  n)  iterations 

to  reduce  the  initial  error  by  a factor  n p > 0.  The  number  of 

multiplications  required  to  reduce  the  error  by  the  same  factor  are 
3 A 

0(n  log  n)  in  two  and  0(n  log  n)  in  three  dimensions. 

In  the  next  section  we  present  results  of  numerical  experiments 
with  MCR  for  the  two-and-three-dimensional  model  problems  for  various 


i 

J 


values  of  a and  h 
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The  storage  and  work  requirements  for  these  problems  are  given  In 
Table  1.  For  the  model  problems,  the  diagonal  elements  of  A are  all 
equal  and  all  the  non-zero  off-diagonal  elements  are  Just  -I's  so  that 
no  storage  is  required  for  A,  and  computing  a matrix-vector  product 
requires  only  N multiplications.  Table  1 also  gives  the  storage  and 
work  requirements  for  a general  problem.  A general  problem  is  one  that 
arises  from  the  application  of  the  simplest  (l.e.,  five-point  in  two  and 
seven-point  in  three  dimensions)  finite  difference  approximations  to  any 
second-order  self-adjoint  elliptic  partial  differential  equation  where  A 
has  the  same  non-zero  structure  as  for  the  corresponding  model  problem. 

A is  still  symmetric,  but  the  off-diagonal  non-zero  entries  need  not  be 
-I's  so  that  the  storage  and  work  costs  are  higher  than  for  the 
corresponding  model  problem. 
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6.  Numerical  Results 

Numerical  experiments  were  performed  to  demonstrate  the  performance 
of  the  MCR  method  on  the  two-and  three-dimensional  model  problems  for 
various  o and  h.  The  solution  was  chosen  to  be 


w(x,y)  - 3e*e^(x-x^){y-y^) 


in  two  and 


w(x,y,z)  - 3e’‘e^e^(x-x^)(y-y^)(z-z^) 


in  three  dimensions.  If  . - 0,  all  the  eigenvalues  of  A are  positive. 
The  first  few  eigenvalues  of  the  differential  operator  -A  in  two 
dimensions  are  19.7,  49.4,  79.0,  98.7,...,  cf.  [8],  so  that  for  o “ 30 
and  90  A has  one  and  three  negative  eigenvalues  respectively. 
Similarly,  in  three  dimensions,  the  eigenvalues  of  -u  are  29.6,  59.2, 
88.8,  108.6,  ...  so  that  for  o • 50  and  1 00  A again  has  one  and  three 
negative  eigenvalues  respectively. 


The  initial  guess  x^  was  taken  to  be  the  zero  vector  and  the  error 


at  Iteration  1 was  defined  as 


error 


f - AXj  II  2/  II  f - AXj^  II  2 


1 /2  -4 

and  was  computed  as  ( (r^,r^) /(r^.r^) ) . We  chose  e « 10  ; 

the  computations  were  carried  out  on  a PDP-IO  (word-length  36  bits)  in 

single  precision. 

Figures  la  and  lb  show  the  error  versus  number  of  iterations  for 
the  model  problems.  Tables  2 and  3 give  the  total  work  required  to 
reduce  the  error  to  10  In  these  experiments  it  turned  out  that 
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MCR  always  chose  the  shorter  Iteration. 

Since  the  difference  approximations  are  only  second-order  accurate, 

we  also  computed  the  number  of  Iterations  required  to  reduce  the  error 
2 

by  a factor  of  1/n  for  various  h and  ..  We  see  from  Figures  2a  and 
2b  that  for  a fixed  problem  (i.e.,  fixed  o)  the  plot  of  the  number  of 
iterations  against  n log  n Is  a straight  line.  This  illustrates  the 
conclusions  of  Theorem  5.1. 

li 
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2 

Two  dimensions  (N-n  ) 

3 

Three  dimensions  (N“n  ) 

Storage 

No.  of  mults./ 

Storage 

No.  of  mults./ 

reqd. 

iteration 

reqd. 

iteration 

Problem 

|aj>c  |a.|<e 

1 a 1 >c  1 a J «c 

1 i - 

1 1 - 

Model 

7N 

8N+2 

ION +4 

7N 

8N+2 

ION +4 

General 

10N-2n 

12N-4n+2 

1 4N  -4  n+4 

llN-3n^ 

14N-6n^+2 

1 6N-6n^+4 

Table  1;  StoraRe  requirements  and  multiplication  counts  for 
N X N systems. 


Mesh- 

No.  of  mults./ 

a 

- 30 

0 

- 90 

width 

iteration 

No.  of 

Total  no. 

No.  of 

Total  no 

h 

la  |>e  la  |<c 

iters. 

of  mults. 

iters. 

of  mults 

1/8 

394 

494 

21 

8,175 

29 

11,327 

1/16 

1,802 

2,254 

52 

93,253 

63 

113,075 

1/32 

7.690 

9,614 

108 

828,597 

131 

1,005,467 

Table  2:  Work  required  to  reduce  the  error  to  10  ^ for 

two-dimensional  model  problems  on  n x n mesh  (h  - l/(n  +1)). 


Mesh- 

width 

h 

No.  of  mults./ 
iteration 

0 

No.  of 
iters. 

- 50 

Total  no. 
of  mults. 

a 

No.  of 
iters. 

- 1 00  1 
Total  no. 
of  mults. 

1/4 

218 

274 

9 

1,907 

8 

1.689  ' 

1/8 

2,746 

3,434 

32 

87, 185 

52 

142,105  1 

1/16 

27,002 

33,754 

71 

1,910,391 

93 

2,504,435 

Table  3:  Work  required  to  reduce  the  error  to  10  ^ for 

three-dimensional  model  problems  on  n x n x n mesh  (h“l/(n+l)). 


mens ions 


Error  versus  number  of  iterations  for  two  and  three-dimensional  model  problems 
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